R Notebook

This is an R Markdown Notebook. When you execute code within the notebook, the results appear beneath the code.

Try executing this chunk by clicking the Run button within the chunk or by placing your cursor inside it and pressing Cmd+Shift+Enter.

library(ape)
library(TreeTools)
library(phangorn)
library(randomcoloR)
library(ggplot2)
library(dplyr)
library(uniqtag)
library(phytools)
library(ggstance)
library(ggtree)
library(grid)
library(gtable)
library(reshape2)
library(ggnewscale)
library(htmlwidgets)

Download global overview tab

raw.data<-read.csv("../../misc/data/GAGA_SPEC_TUBE_ECOL-GlobalOverview-210322.tsv",sep="\t",strip.white = T,stringsAsFactors = F,dec = ",",na.strings = c("","#DIV/0!"),skip = 2)

Load from google

# subset to entries with a Species name from the Assembly Overview
data<-subset(GlobalOverview,!is.na(Species..from.Assembly.Overview.) & !is.na(ENTRY.FROM.ASSEMBLY.OVERVIEW))

# rename Myrmoxenus to Temnothorax as the former is no proper genus
data$Species..from.Assembly.Overview.<-gsub("Myrmoxenus","Temnothorax",data$Species..from.Assembly.Overview.)

# Extract genus information 
data$genus<-gsub(" .*","",data$Species..from.Assembly.Overview.)

# Rename duplicate genera to Genus-sp1, Genus-sp2, etc.
data$taxID<-make_unique(data$genus,sep="-sp")

# show which genera are duplicated
table(gsub("-.*","",data$genus[duplicated(gsub(" .*","",data$Species..from.Assembly.Overview.))]))+1

   Acromyrmex Aphaenogaster          Atta    Camponotus      Carebara Crematogaster      Diacamma  Ectomomyrmex       Formica  Harpegnathos        Lasius 
            5             3             2             6             7             3             2             2             6             2             8 
   Leptogenys        Messor      Myrmecia       Myrmica  Odontomachus      Pheidole   Proceratium    Proformica   Rossomyrmex    Solenopsis      Tapinoma 
            3             2             3             3             3             5             2             2             2             2             2 
  Temnothorax   Tetramorium  Trachymyrmex 
            8             5             3 

# Downloaded ant tree from https://www.antwiki.org/wiki/Phylogeny_of_Formicidae
## The tree contains a few typos that I corrected manually (e.g. Paratrachina instead of Paratrechina, Myrmecosystus instead of Myrmecocystus)
treePaste<-paste(scan("../../misc/data/antTree.tre", what="character", sep=" "),collapse=" ")

Read 1 item

# expand the nodes for those species that we have multiple times
expansions<-table(data$genus[duplicated(data$genus)])+1
for (i in 1:length(expansions)){
  treePaste<-gsub(paste("(",names(expansions)[i],":)",sep=""),paste("(",paste(rep("\\1",expansions[i]),collapse=":0,",sep=""),":0):",sep=""),treePaste)
}

# load manipulated tree as tree object
raw.tree<-read.tree(text=treePaste)
raw.tree$tip.label<-make_unique(raw.tree$tip.label,sep="-sp")

# show labels for figuring out subfamilies
plot(raw.tree,show.tip.label = F)

92 branch length(s) NA(s): branch lengths ignored in the plot

nodelabels(frame =  "none")

# subset to extant ants, Martialis has to be node 1
tree<-Subtree(Preorder(raw.tree), 427)
grep("Martialis",tree$tip.label)

[1] 1

tree$tip.label<-gsub("some(.*[2-9])|some(.*)1","\\1\\2",tree$tip.label)

# plot with subfamilies

p<-ggtree(tree,size=.1)+xlim(0,40)+geom_tiplab(align=T,size=1)+ geom_text(aes(label=node),size=1,col="darkred")

Duplicated aesthetics after name standardisation: size

# Download species counts for all genera and other info from antwiki
#curl https://www.antwiki.org/wiki/images/a/ad/AntWiki_Valid_Genera.txt | LC_ALL=C sed 's/[^[:blank:][:print:]]//g'  > ../misc/data/AntWiki.clean.tsv

speciescount<-read.csv("../../misc/data/AntWiki.clean.tsv",sep="\t")
speciescount$Tribe[speciescount$Tribe==""]<-paste("tr. ",speciescount$Subfamily[speciescount$Tribe==""],sep="")

dataMerge<-merge(data,speciescount,by.x="genus",by.y="TaxonName",all.x=T,all.y=T)

dataMerge <- dataMerge %>%
  select(taxID, everything())

# give name to all unsequenced genera
dataMerge$taxID[is.na(dataMerge$taxID)]<-dataMerge$genus[is.na(dataMerge$taxID)]
# set species Count to NA for all multiple genera species
dataMerge$SpeciesCount[grep("-sp[2-9]",dataMerge$taxID,value = F)]<-NA

add data to tree

p<-ggtree(tree,size=.2,layout = "circular")+xlim(0,40)
#p$data$label<-make_unique(p$data$label,sep = "-sp")
#p$data$label[grep("NA.",p$data$label,ignore.case = F)]<-NA
dataMerge$DisplayName<-gsub(" \\(.*","",dataMerge$Species..from.Assembly.Overview.)
p2<-p %<+% dataMerge 

antGenera.tree <- p2 +
          geom_tiplab(align=T,size=1.2,offset=.5,linesize=0)+
          geom_tippoint(aes(size = SpeciesCount,fill=Tribe),pch=21,alpha=.5) + 
          guides(fill="none",color="none")+
          theme_tree()+
          theme(legend.position = c(.9,.85),legend.background = element_blank())

Duplicated aesthetics after name standardisation: sizeDuplicated aesthetics after name standardisation: size

antGenera.tree.w.GAGA <- p2 +
          geom_tiplab(aes(label=DisplayName,col=Tribe,alpha=ifelse(is.na(DisplayName),0,1)),align=T,size=1.8,offset = 8,linesize = 0)+
          geom_tiplab(align=T,size=1.2,offset=.5,linesize=0)+
          geom_tippoint(aes(size = SpeciesCount,fill=Tribe),pch=21,alpha=.5) + 
          guides(fill="none",color="none",alpha="none")+
          theme_tree()+
          theme(legend.position = c(.9,.85),legend.background = element_blank())

Duplicated aesthetics after name standardisation: sizeDuplicated aesthetics after name standardisation: sizeDuplicated aesthetics after name standardisation: sizeDuplicated aesthetics after name standardisation: size

antGenera.tree
ggsave("../../misc/plots/antTree.pdf",height=10)

Saving 7.29 x 10 in image

antGenera.tree.w.GAGA
ggsave("../../misc/plots/antTree.GAGA.pdf",height=10)

Saving 7.29 x 10 in image

# How many GAGA species are shown?
sum(!is.na(dataMerge$DisplayName))

[1] 163

GAGA.tree<-keep.tip(tree,antGenera.tree.w.GAGA$data$node[!is.na(antGenera.tree.w.GAGA$data$DisplayName)])

GAGA.tree.plot<-ggtree(GAGA.tree,size=.2) %<+% dataMerge 
GAGA.tree.plot2<-GAGA.tree.plot+geom_tiplab(aes(label=DisplayName,col=ifelse(grepl("GAGA",ENTRY.FROM.ASSEMBLY.OVERVIEW),"GAGA","other")),size=1.5,show.legend=F)+scale_color_manual(name="",values = c("black","grey60"))

GAGA.tree.plot2

GAGA.tree.plot3

# extract data from tree for facet plotting

facetData<-as.data.frame(GAGA.tree.plot2$data[GAGA.tree.plot2$data$isTip==T,])
facetData$id<-facetData$label
facetData$Genome.size<-as.numeric(facetData$Genome.size)
facetData$Scaffold.N50<-as.numeric(facetData$Scaffold.N50)
facetData$Contig.N50<-as.numeric(facetData$Contig.N50)
facetData$Number.of.scaffolds<-as.numeric(facetData$Number.of.scaffolds)

# extract BUSCO scores
facetData$BUSCO.hym.C<-as.numeric(gsub("C:(.*?)%.*","\\1",facetData$BUSCO.Hymenoptera))
facetData$BUSCO.hym.F<-as.numeric(gsub(".*F:(.*?)%.*","\\1",facetData$BUSCO.Hymenoptera))
facetData$BUSCO.hym.M<-as.numeric(gsub(".*M:(.*?)%.*","\\1",facetData$BUSCO.Hymenoptera))

facetData$BUSCO.euk.C<-as.numeric(gsub("C:(.*?)%.*","\\1",facetData$BUSCO.Eukaryota))
facetData$BUSCO.euk.F<-as.numeric(gsub(".*F:(.*?)%.*","\\1",facetData$BUSCO.Eukaryota))
facetData$BUSCO.euk.M<-as.numeric(gsub(".*M:(.*?)%.*","\\1",facetData$BUSCO.Eukaryota))


dd <- data.frame(id=facetData$label, GS=round(facetData$Genome.size/1000000),N50=facetData$Scaffold.N50,nscf=facetData$Number.of.scaffolds,B.hym.C=facetData$BUSCO.hym.C,B.hym.F=facetData$BUSCO.hym.F,B.hym.M=facetData$BUSCO.hym.M,B.euk.C=facetData$BUSCO.euk.C,B.euk.F=facetData$BUSCO.euk.F,B.euk.M=facetData$BUSCO.euk.M,tribe=facetData$Tribe,PacBio=facetData$PacBio.1,stLFR=facetData$stLFR.2,HiC=facetData$Hi.C.1,RNAseq=facetData$RNA..additional.samples.still.under.seq..,Finalized=facetData$Current.state,wgs=facetData$Short.reads..WGS..1,Gid=facetData$GAGA.ID,status=facetData$Current.state)
dd$tribeColor <- cladeColor[match(dd$tribe, names(cladeColor))]

dd$PacBio[!is.na(dd$PacBio)]<-1
dd$PacBio[is.na(dd$PacBio)]<-0
dd$stLFR[!is.na(dd$stLFR)]<-1
dd$stLFR[is.na(dd$stLFR)]<-0
dd$HiC[!is.na(dd$HiC)]<-1
dd$HiC[is.na(dd$HiC)]<-0
dd$wgs[!is.na(dd$wgs)]<-1
dd$wgs[is.na(dd$wgs)]<-0

rownames(dd)<-dd$id

dd$RNAseq[!is.na(dd$RNAseq)]<-3
dd$RNAseq[is.na(dd$RNAseq)]<-0

dd$Finalized[grep("Final",dd$Finalized)]<-2
dd$Finalized[dd$Finalized!=2]<-0

dd$DisplayName<-facetData$DisplayName

dd$node<-facetData$node
dd$tech<-as.factor(paste(dd$PacBio,dd$HiC,dd$stLFR,sep=""))
levels(dd$tech)[levels(dd$tech)=="000"]<-"non-GAGA"
levels(dd$tech)[levels(dd$tech)=="100"]<-"PacBio"
levels(dd$tech)[levels(dd$tech)=="001"]<-"stLFR"
levels(dd$tech)[levels(dd$tech)=="101"]<-"P/s"
levels(dd$tech)[levels(dd$tech)=="110"]<-"P/H"
levels(dd$tech)[levels(dd$tech)=="111"]<-"P/H/s"


ddHeatmap<-dd[c("id","tribeColor","tribe","PacBio","stLFR","HiC","wgs","Finalized","RNAseq")]

#GAGA.tree.plot4<-facet_plot(GAGA.tree.plot3, 'Genome Size', data = dd, geom=geom_barh, mapping = aes(x = GS,fill=tribeColor), 
 #               stat='identity' )+ theme_tree2() + scale_fill_identity()

GAGA.tree.plot.heatmap<-gheatmap(GAGA.tree.plot3+ylim(0,165), ddHeatmap[,c(4:9)], offset=7, colnames=T, legend_title="Tech",  width = .25,colnames_position = "top", colnames_offset_y = 1,colnames_angle=65,hjust=0,font.size = 2.8)+
                  scale_fill_manual(values = c("grey90","lightgreen","darkred","lightblue"))+
                  theme(legend.position = "none")

Scale for 'y' is already present. Adding another scale for 'y', which will replace the existing scale.
Scale for 'fill' is already present. Adding another scale for 'fill', which will replace the existing scale.

GAGA.tree.plot.heatmap
ggsave("../../misc/plots/GAGA.tree.technique.pdf",GAGA.tree.plot.heatmap,width=6,height=10)

ddHeatmap2<-dd[c("id","tribeColor","tribe","B.hym.C","B.euk.C")]
colnames(ddHeatmap2)[4:5]<-c("Hymenopt.","Eukaryota")
row.names(ddHeatmap2)<-ddHeatmap2$id

heatcol <- c("steelblue","cyan", "red")
GAGA.tree.plot.busco<-gheatmap(GAGA.tree.plot3+ylim(0,165), ddHeatmap2[,c(4:5)], offset=7, colnames=T, legend_title="Tech",  width = .215,colnames_position = "top", colnames_offset_y = 1,colnames_angle=65,hjust=0,font.size = 2.5)+scale_fill_gradientn(colours = rev(heatcol),na.value = "grey90",name="% comp. Busco")+theme(legend.position = c(.2,.8),legend.background = element_blank())

Scale for 'y' is already present. Adding another scale for 'y', which will replace the existing scale.
Scale for 'fill' is already present. Adding another scale for 'fill', which will replace the existing scale.

GAGA.tree.plot.busco
ggsave("../../misc/plots/GAGA.tree.busco.pdf",GAGA.tree.plot.busco,width=6,height=10)

GAGA.tree.plot.N50<-facet_plot(GAGA.tree.plot.heatmap+new_scale("color"), 'scfN50 (Mb)', data = dd, geom=geom_point, mapping = aes(x = N50/1000000,color=tech),size=2)+ 
  theme_tree2()+
  scale_color_brewer(palette = "Set2" ,guide = guide_legend(nrow=2,label.position = "top"))+
  theme(legend.position = c(.9,.9),
        legend.background = element_rect(color="black",fill=rgb(0,0,0,0)),
        legend.margin = margin(0,2,2,2,"mm"),
        legend.spacing = unit(0,"mm"),
        legend.title = element_blank(),
        legend.key.height = unit(0,"mm"),
        legend.key.width = unit(0,"mm"),
        legend.box.spacing = unit(0,"mm"))+
        guides(fill="none")+
  ylim(0,170)

Scale for 'y' is already present. Adding another scale for 'y', which will replace the existing scale.

GAGA.tree.plot.N50

# change rel_width of panels
## see https://groups.google.com/g/bioc-ggtree/c/tZ0qkluBeGU/m/nkfWC9ixCwAJ

gt = ggplot_gtable(ggplot_build(GAGA.tree.plot.N50))

Removed 26 rows containing missing values (geom_point).

#gtable_show_layout(gt) # will show you the layout - very handy function
#gt # see plot layout in table format
gt$layout$l[grep('panel-1', gt$layout$name)] # you want to find the column specific to panel-2

[1] 5 7

gt$widths[5] = 1.9*gt$widths[5] # in this case it was colmun 7 - reduce the width by a half
#grid.draw(gt) # plot with grid draw

grid.draw(gt)
ggsave("../../misc/plots/GAGA.tree.n50.pdf",grid.draw(gt),width=10,height=9)

ggsave("../../misc/plots/GAGA.assemblyQuality.pdf",assembly.quality.plot2,width=6,height=6)

ggsave("../../misc/plots/GAGA.assemblyQuality.boxplots.pdf",GAGA.quality.boxplots,width=6,height=6)

dd$Subfamily<-facetData$Subfamily.y
dd$DisplayNameU<-make_unique(dd$DisplayName,sep="")
dd$DisplayNameU<- factor(dd$DisplayNameU, levels = dd$DisplayNameU[rev(order(dd$Subfamily,dd$GS))])